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Abstract 

We report an empirical determination of the probability density functions -Pdata( r ) of the number r of 
earthquakes in finite space-time windows for the California catalog, over fixed spatial boxes 5x5 km 2 
and time intervals dt = 1, 10, 100 and 1000 days. We find a stable power law tail -Pdata( r ) ~ l/r 1+ ^ 
with exponent \i ~ 1.6 for all time intervals. These observations are explained by a simple stochastic 
branching process previously studied by many authors, the ETAS (epidemic-type aftershock sequence) 
model which assumes that each earthquake can trigger other earthquakes ("aftershocks"). An aftershock 
sequence results in this model from the cascade of aftershocks of each past earthquake. We develop 
the full theory in terms of generating functions for describing the space-time organization of earthquake 
sequences and develop several approximations to solve the equations. The calibration of the theory to 
the empirical observations shows that it is essential to augment the ETAS model by taking account of 
the pre-existing frozen heterogeneity of spontaneous earthquake sources. This seems natural in view of 
the complex multi-scale nature of fault networks, on which earthquakes nucleate. Our extended theory is 
able to account for the empirical observation satisfactorily. In particular, the adjustable parameters are 
determined by fitting the largest time window dt = 1000 days and are then used as frozen in the formulas 
for other time scales, with very good agreement with the empirical data. 



*Electronic address: sornette@moho.ess.ucla.edu 
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I. INTRODUCTION 



Many papers purport to characterize the space-time organization of seismicity in different regions 
of the world. Recent claims of universal laws for the distribution of waiting times and seismic rates 
between earthquakes have derived from the analyses of space-time windows The flurry 

of interest from physicists comes from their fascination with the self-similar properties exhibited 
by seismicity (Gutenberg-Richter power law of earthquake seismic moments, Omori decay law 
of aftershock rates, fractal and multifractal space-time organization of earthquakes and faults) 
logether with the development of novel concepts and techniques that may provide new insights 

The interest is no less vivid among seismologists and geophysicists in characterizing the space- 
time properties of seismicity, because it allows them to understand the dynamics of plate motion 
(at large scales), to constrain the interaction between faults, as well as to develop better hazard as- 
sessment. Recently, an additional incentive is provided by the development of forecasting models of 
seismicity, for instance within the RELM (Regional Earthquake Likelihood Models: www.relm.org) 
project in Southern California. In the RELM project, a forecast is expressed as a vector of earth- 
qU a k e rates specaed for each mU M- dlm e nslonal bi n g where a bin is defined by an interval 
of location, time, magnitude and focal mechanism and the resolution of a model corresponds to 
the bin sizes. Then, expectations and likelihoods can be estimated and used for the comparison 
between different forecasting methods. 

A fundamental issue is testing models' prediction is to take into account so-called aftershock 
clustering. In one way or another, most if not all models use some form of declustering approach 
to remove the effect of aftershocks which otherwise dominate and obscure the desired information 
about the model's performance [8|. Then, with such declustered catalogs, the likelihood of forecasts 
are estimated using Poissonian probabilities. But, if the catalog is only partially declustered (which 
it will most probably be as there are no agreed upon fully efficient method of declustering), then 
our contribution in this paper is to show that the distribution of event numbers should present 
a tail much more heavy than predicted by the Poissonian statistics and to propose a theoretical 
explanation for it. Pisarenko and Golubeva [l| introduced a model to decluster catalogs by so-to- 
say Poisson "with random parameter," which resulted in a law with slowly decreasing probabilities 
(actually a stable Levy law) for the distribution of rates. 
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We improve on preceding results on several points. Our first contribution is to show that the 
heavy tail nature of the distribution of seismic rate is intrinsic to a class of generic models of 
triggered seismicity. Specifically, our theory is based on a simple model of earthquake triggering, 
in which future seismicity is a conditional Poisson process, with average rates (or Poisson intensity) 
conditioned on past seismicity. We show that the exponential Poisson rate is renormalized into a 
power law tail by the mechanism involving a cascade of earthquake triggering. Our theory thus 
provides a prediction for the distribution of seismic rates in space-time bins in the form of a power 
law tail distribution. Our second contribution is to show that our prediction is verified by empirical 
seismic rates in Southern California over more than two decades. In addition, our theory accounts 
well for the evolution of the distribution of seismic rates as a function of the time window size from 
1 day to 1000 days. 

This implies that spontaneous fluctuations of the number of triggered earthquakes in space-time 
bins may be simply due to the cascades of triggering processes, which lead to dramatic departures 
from the Poisson model used as one of the building block of standard testing procedures. Account- 
ing for the intrinsic heavy tail nature of the distribution of seismic rate may explain, we believe, 
many of the contradictions and rejections of models assessed on the basis of Poisson statistics of 
so-called declustered catalogs. This also suggests the need for fundamentally different earthquake 
prediction models and testing methods. Our results also offer a simple alternative explanation to 
so-called universal laws in terms of cascades of triggered earthquakes: our proposed framework 
explains the observed power law distributions of seismic rates from the fundamentals of seismicity 
characterized by a few exponents. 



II. THE EPIDEMIC-TYPE AFTERSHOCK SEQUENCE (ETAS) BRANCHING 
MODEL OF EARTHQUAKES WITH LONG MEMORY 

We study the general branching process, called the Epidemic- Type Aftershock Sequence (ETAS) 
model of triggered seismicity, introduced by Ogata in the present form and by Kagan and 
Knopoff in a slightly different form and whose main statistical properties are reviewed in 
For completeness and in order to fix notations, we recall its definition and ingredients used in our 
analysis that follows. In this model, all earthquakes are treated on the same footing and there 
is no distinction between foreshocks, mainshocks and aftershocks, other than from retrospective 
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human-made classification. The advantage of the ETAS model is its conceptual simplicity based 
on three independent well-found empirical laws and its power of explanation of other empirical 
observations (see for instance and references therein). 



14. 



3, 



and has in addition 



The ETAS model belongs to a general class of branching processes 
the property that the variance of the number of earthquake progenies triggered in direct lineage from 
a given mother earthquake is mathematically infinite. Moreover, a long-time (power law) memory 
of the impact of a mother on her first-generation daughters describes the empirical Omori law for 
aftershocks. These two ingredient s tog ether with the mechanism of cascades of branching have been 
shown to give rise to subdiffusion [ly, Ujj and to non mean-field behavior in the distribution of the 
total number of aftershocks per mainshock, in the distribution of the total number of generations 
before extinctions lal and in the distribution of the total duration of an aftershock sequence before 

n 

extinction |19l |. 

In the ETAS model, each earthquake is a potential progenitor or mother, characterized by its 
conditional average number 

N m = Kn{m) (1) 
of children (triggered events or aftershocks of first generation), where 

H{m) = 10 a(m " mo) , (2) 

is a mark associated with an earthquake of magnitude m ^ mo (in the language of "marked point 
processes"), k is a constant factor and mo is the minimum magnitude of earthquakes capable 
of triggering other earthquakes. The meaning of the term "conditional average" for N m is the 
following: for a given earthquake of magnitude m and therefore of mark fi(m), the number r of its 
daughters of first generation are drawn at random according to the Poissonian statistics 

^) = f^ = ^e-. (3) 

N m is the expectation of the number of daughters of first generation, conditioned on a fixed 
magnitude m and mark /i(m). The expression (J2J) for //(m) is chosen in such a way that it 
reproduces the empirical dependence of the average number of aftershocks triggered directly by 
an earthquake of magnitude m (see j3| and references therein). Expression (0) with (j2J) gives the 
so-called productivity law of a given mother as a function of its magnitude. The challenge of our 
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present analysis is to understand how the exponential distribution (jSJ) is changed by taking into 
account all earthquake triggering paths simultaneously and at all possible generations. 

The ETAS model is complemented by the Gutenberg- Richter (GR) density distribution of earth- 
quake magnitudes 

p(m) = b ln(10) l()- b ( m - m °) , m ^ m , (4) 

such that p(x)dx gives the probability that an earthquake has a magnitude equal to or larger 
than m. This magnitude distribution p(m) is assumed to be independent of the ma gnit ude of the 
triggering earthquake, i.e., a large earthquake can be triggered by a smaller one [r! I2 fl| . 

Combining (j3J) and (J2J) shows that the earthquake marks fi and therefore the conditional average 
number N m of daughters of first generation are distributed according to the normalized power law 

vM = -y+j, l</i<+oo, y = b/a. (5) 



For earthquakes, b « 1 and 0.5 < a < 1 giving 1 < 7 < 2 (see 2l( for a review of values quoted 
in the literature and their implications). This range 1 < 7 < 2 implies that the mathematical 
expectation of ji and therefore of N m (performed over all possible magnitudes) is finite but its 
variance is infinite (the marginal case a = 1 leading to 7 = 1 requires the existence of an upper 
magnitude cut-off j2l|). 

For a fixed 7, the coefficient k then controls the value of the average number n (or branching 
ratio) of children of first generation per mother: 

7 

n = (N m ) = k(h) = k , (6) 

7-1 

where the average (N m ) is taken over all mothers' magnitudes drawn from the GR law. Recall 
that the values n < 1, n = 1 and n > 1 correspond respectively to the sub-critical, critical and 
super-critical regimes. 

The next ingredient of the ETAS model consists in the specification of the space-time rate 
function N m $(r — Ti,t — ti) giving the average rate of first generation daughters at time t and 
position r created by a mother of magnitude m ^ itlq occurring at time ti and position r^. 

<S>(x,t) = <S>(t)<f>(x). (7) 



The time propagator $(t) has the Omori law form 

6d 
(c + t) 



= TTT^ITe Hi*) ( 8 ) 



where H(t) is the Heaviside function, < 9 < 1, c is a regularizing time scale that ensures that 
the seismicity rate remains finite close to the mainshock. The time decay rate (JBJ) is called the 



"direct Omori law" 



12J, |22j. Due to the process of cascades of triggering by which a mother 



triggers daughters which then trigger their own daughters and so on, the direct Omori law (jSJ) is 



which is the one observed 



renormalized into a "dressed" or "renormalized" Omori law 
empirically. The analysis below will retrieve and extend this result. 
The space propagator is given by 

, . r] d v 

0(:E)= 27T(x 2 + rf 2 )(^ 2 )/ 2 ' ( ) 

For our comparison with the empirical data, we shall consider the epicenter position of earthquakes, 
that is, the 2D-projection on the earth surface of the real 3D-distribution of earthquake hypocenters. 
Numerical implementations of the theory developed below will thus be done in 2D but it is easy 
to generalize to 3D if/when the empirical data will be of sufficient quality to warrant it. 
In the following, we will make use of the Laplace transform of the Omori law 



$(u)=/ $(t)e- ut dt = 6(cu) e e cu T(-6,cu) (10) 
Jo 

and of its asymptotic behavior 

$ _1 (u) ~ 1 + T(1 - 9)(cu) e , cu<l. (11) 
The Fourier transform of the space propagator @ will also be useful: 



In particular 



e- dq fo = l), i{q) = e- d *{l + dq) = 3) . (13) 



The last ingredient of the ETAS model is to assume that plate tectonic motion induces spon- 
taneous mother earthquakes, which are not triggered by previous earthquakes, according to a 
Poissonian point process, such that the average number of spontaneous mother earthquakes per 
unit time and per unit surface is g. In the ETAS branching model, each such spontaneous mother 
earthquake then triggers independently its own space-time aftershocks branching process. 
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It is a well-established facts that, at large scale, earthquakes are preferentially clustered near 
the plate boundaries while, at smaller scales, earthquakes are found mostly along faults and close 
to nodes between several faults (^j). It is natural to extend the ETAS model to allow for the 
heterogeneity of the spontaneous earthquake sources g reflecting the influence of pre-existing fault 
structures, some rheological heterogeneity and complex spatial stress distributions. We get some 
guidelines from the distribution of the stress field in heterogeneous media and due to earthquakes 



24 1 which should be close to a Cauchy distribution or probably more generally to a power law 
25, Q (see also Chap. 17 of H) 

The simplest prescription is thus to assume that g is itself random and distributed according to 

A' to)- 

where (g) is then statistical average of the random space-time Poissonian source intensity g. In 
the numerical applications below, we shall use the form 

f 5 (x) = °-— 1 + | (5>0), (15) 



5 V 5 

The value 5 = gives the same tail as the Cauchy distribution advocated in [24] for the stress field. 
We have considered other functions, such as half-Gaussian, exponential, half-Cauchy but none of 
them give satisfactory fits to the data (see below). The parametrization (fTB]) with 5 > allows us 
to have only a single scale (g) controlling the typical fluctuation of the random sources. We have 
found that only slightly positive values of 5 (corresponding to tails a little fatter than the Cauchy 
law) gives reasonable fits to the data (see below). It is interesting to observe that the data on the 
distribution of seismic rates thus seems to constrain significantly the fractal distribution of seismic 
sources. 



III. GENERATING PROBABILITY FUNCTION (GPF) OF EARTHQUAKES 
BRANCHING PROCESS 

In this section, we describe the statistical properties of earthquake branching processes using 
the technology of generating functions. 

First, let us recall the GPF of the total number Ri of the first-generation aftershocks of a mother 



8 



event of magnitude m can be easily obtained as 

e*«(*-i) ; (1 6 ) 

using the fact that the rate of first-generation aftershocks is Poissonian according to (J3J). In this 
expression (fTBj). n and fi are given by their definition in (JTJ) and ©• 

Averaging over the random parameter /i gives the GPF of the number Ri of first generation 
aftershocks triggered by a mother aftershock of arbitrary magnitude 

G(z) =7K 7 (l-z) 7 r(- 7 ,«(l-^)). (17) 

Note that is nothing but the branching ratio n defined above in equation (jUJ). Knowing the 
expression (|T7j) for the GPF G(z), one finds that the corresponding probabilities of the random 
numbers Ri are equal to 

Pi(r) =Pr{ J R 1 = r} = 7— r(r - 7, «) , (18) 

T! 

which have the following asymptotics 

Pi(r) - ^ = n 7 7^(7- l)^- 7 - 1 (r>l). (19) 

Expression (JT9)l implies that, for 1 < 7 < 2, the variance of the random number R\ is infinite. For 
7 > 2, the variance is finite and is equal to 

Figure 1 shows the probabilities (jTHj) and their power law asymptotics for n = 1, 7 = 1.25 and 
7 = 3. 

Let us now consider the set of independent space-time aftershock branching processes, triggered 
by spontaneously arising mother earthquakes. Due to the independence between each sequence 
triggered by each spontaneous event, it is easy to show that the GPF of the number of events 
(including mother earthquakes and all their aftershocks of all generations), falling into the space- 
time window {[t, t + r] x S} is equal to 

e sp (z,r,S)=e-e L ^ (21) 
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where 



oo 

roo 



/•oo p r 

L{z, t,S)= dt dx[l - Q{z, t, t, S; x)}+ 

— oo 

oo 

J T dt J J dx[l-Q(z,t,S;x)][l - I s (x)]+ (22) 

— oo 

J T dt J J dx[l - zQ(z, t, S; x)] . 
s 

The three above summands have the following transparent geometrical meaning. 

• The first summand describes the contribution to the GPF sp from aftershocks triggered 
by mother earthquakes that occurred before the time window (i.e. at instants if such that 
t' < t) (positions 1 in Fig. 2). The corresponding GPF Q(z,t — t',r,S;x) of the number of 
aftershocks triggered inside the space-time window {[t, t + r] xS} by some mother event that 
occurred at time t' satisfies the relation 

B(z, t, r, S; x) = G[l - V{z, t, t, S; x)} , (23) 

where the auxiliary function ty(z, t, r, S; x), describing the space-time dissemination of after- 
shocks triggering by some mother event, is equal to 

V(z,t,r,S\ x) = 

$(sr, t)®[l- Q(z, t, t, S; x)} + $(x, t + r)® [I- Q{z, r, S; x)}+ (24) 
(1 - z)$(oj, t + r)® I s {x)Q(z, r, S; x) . 

$(a? — x',t'), which has been defined in flj|), is the probability density function (pdf) of the 
position x' and instant if of some first generation aftershock, triggered by the mother event, 
arising at the instant t = and at the point x. The function Is{ x ) i n f)24|) is the indicator 
of the space window S and G(z) in (}23|) is the GPF of the number R\ of first generation 
aftershocks, triggered by some mother earthquake. G(z) given in (j!7)) is common to all 
mother earthquakes and to all aftershocks. 

• The last two terms in expression ()22j) (positions 2 and 3 in Fig. 2) describe the contribution 
of aftershocks triggered by earthquakes, occurring inside the time window (i.e., t' G [t, t + r]). 
The second term (position 2 in Fig. 2) corresponds to the subset spatially outside the domain 
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S. The third term (position 3 in Fig. 2) corresponds to the subset spatially inside the domain 
S. These last two terms in expression (|22|) depend on the GPF 



Q(z,r,S;x) = S(z,t = 0,r,S;x) (25) 

of the numbers of aftershocks triggered till time r inside the space window S by some mother 
event arising at the instant t — and at the point x. It follows from (J23j) and (J24j) that it 
satisfies the relations 

Q(z,r,S;x) = G[l - $(z,T,S;x)} (26) 

and 

ty(z,r,S: x) = 

(27) 

$(x, t) ® [1 - Q(z, t, S; x)] + (1 - z)$(x, r) <g> J 5 (a;)e(z, r, 5; aj) . 
In addition, we shall need the GPF 

G(z,S;x) = 6(2, r = oo,S;x) (28) 

of the total numbers of aftershocks triggered by some mother earthquake inside the area S. As 
seen from (|26|) and (|27|). it satisfies the relations 

G(z,S;x) = G[l-V(z,S;x)] (29) 

and 

$?(z, S;x) = l- 4>(x) <g> 6(^, 5; sc) + (1 - z)0(a;) <g> 7 5 (aj)e( i z, 5; a:) . (30) 

Taking into account the distribution of the source intensities g amounts to averaging equation 
(J2*T|) over g weighted with the statistics ([HJ). This gives 

e sp (z,T;S) = f[(g)L(z,T,S)}, (31) 

where f(u) is the Laplace transform of the pdf f(x). For the specification (fTHj) . expression (|3"Tf 
becomes 

= (l + 5)(<k0 1+ V u r(-l-5,foO. (32) 
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IV. AVERAGES AND RATES OF AFTERSHOCKS WITHIN THE SPACE-TIME WIN- 
DOW {[t,t + r] x S} 



Before discussing the properties of the distributions of aftershocks, we consider their simplest 
statistical characteristics, namely the averages and rates of different kinds of aftershocks. This 
introduces the relevant characteristic scales in the time and in the space domains, which are found 
inherent to the space-time branching processes. This also suggests the natural "large time window 
approximation" used and tested below within the more general probabilistic treatment. 



A. Average of the total number of events in the space time window {[t, t + r] x S} 

Let us first calculate the average the total number of events inside the space-time window given 



by 



(Rs P (r,S)) 



dQ sp (z,r, S) 



dz 

It follows from (|3Tj) and (j22J) that it is equal to 



(33) 

2 = 1 



(R sp (r,S)) = (R out {r,S)) + (R{t,S)) + (q)St, (34) 

where (-R out (r, 5)) is the average number of aftershocks triggered by spontaneous "mother" earth- 
quake sources that occurred before time t (positions 1 in Fig. 2), (R(r, S)) is the average of number 
of aftershocks triggering by spontaneous earthquake sources that occur within the time interval 
[t, t+r] (positions 2 and 3 in Fig. 2) and (g)Sr is the average of number of spontaneous earthquakes 
inside the space-time window. Here and everywhere in the following, S is the area of the spatial 
domain S and thus St is the space-time volume associated with the space-time window. 
In what follows, it will be useful to introduce the rate of events 

iV !p(T , 5H ^M!>, (35) 

with 

iV sp (r, S) = N out (r, S) + N(r, S) + (g}S . (36) 
Using Eq. (|23 jl . ([21)1 and Eq. l|2fy . (|27jl . one shows that 

77 

N(t,S) = — (£?> S - N ont (r,S) . (37) 
12 



where iV ou t(r, S) satisfies the equation 

AU(t, S) - n N mt (r, S) <g> $(r) = ^^a(r) , (38) 
and n is the branching ratio defined in (JHJ). Here and below, the following notation is used 

a(r) = 1°° Ht) dt = . (39) 

Substituting (|3*7)l into (|3*^jl gives the obvious equality 

iV sp (r,<S) = -^, (40) 

1 — 71 

which implies that, due to the cascade of earthquake triggering processes, the average of the total 
number (i? sp (r, S)) of events is amplified by the factor 1/(1 — n) compared with the average number 
(q)tS of earthquake sources. This factor 1/(1 — n) has a simple intuitive meaning J^J: one event 
gives on average n daughters in direct lineage; each of these first-generation daughters give n grand- 
daughters, the average number of grand-daughters is thus n 2 , and the reasoning continues over all 
generations. Summing over all generations, the total number of events triggered by a given source 
plus the source itself is 1 + n + n 2 + n 3 + which sums to 1/(1 — n). 

B. Impact of mother earthquake sources occurring before the time window 

Here, we use previous and other related relations in the goal of estimating the contribution of 
the different terms in the r.h.s. of Eq. (|22j) to the GPF sp (z, r, S). Our goal is to prepare and 
check for approximations that will be used below. 

For instance, we shall assume that the contribution of the first term in (J22)) . which is responsible 
for aftershocks triggered by earthquakes occurring before time t (i.e., outside the time window 
[t, i + r]), is negligible if the corresponding relative events rate obeys the following condition 

N out (r,S) 

Nout{T) = j- q , < 1 • (41) 
To check when this condition holds, notice that, due to (J3~%j) . fl4*Uj) . A/" ou t(T) is solution of 

Kut(r) - nKut(r) <g> $(r) = na(r) . (42) 
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Applying the Laplace transform to both sides of this equation gives 



n- 



Using the asymptotic formula (fTT|) . we obtain 



Aout(w) 



u[l - n&(u)) 



n(ciu) e 1 



1 + (ciu) 



where 



r(i 



1/0 



C'l 



n 



(43) 



(44) 



(45) 



is a characteristic time-scale of aftershock branching processes separating a l/t 1-6 * law at t < c\ 
from a l/t 1+£l law at t > cj for the decay with time of the average aftershock rate triggered by a 
single mother earthquake [12, 
~ 10.5 hours. 



221 ] . For instance, if c = 2 min, n = 0.9, 6=1/2 then ci ~ 628 min 



Taking the inverse Laplace transform of Eq. 

ACt(r) =nE, 



gives 
r 
ci 



where Eq(x) is the Mittag-Leffler function defined by 



£e(-a;) 



.r 



y° 



dy 



sin7r^ , 

7r jo y ~\~ x -\- zxy cos tiO 



(x > 0, < < I) 



The following asymptotic property holds: 

Ee(-x) 



1 



x rri 



[X 



oo 



(46) 



(47) 



(4* 



In addition, 

Ei/2(-x) = e x2 erfc x . (49) 

Figure 3 plots the exact rate Af out (r), its fractional approximation (jjH)) and corresponding 
asymptotics derived from 



n 



out 



r l- 



(50) 



for n = 0.9 and 9 = 1/2. One can observe that the asymptotic result fl5()D is rather precise even if 
r is close to ci. Eq. (jHilj) means that, if 

r » ci , (51) 

then one can neglect the contribution of aftershocks triggered by the spontaneous earthquake 
sources occurring before the time window [t,t + t\. The remark will be used in our following 
investigation and the condition (J5T| will be refered to as the "large time window approximation." 
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C. Impact of mother earthquake sources occurring inside the space-time window {[t, t + 
r]xS} 



Let calculate the contribution to the rate (|36|) of events corresponding to the last term of (|22j) 
(position 3 in Fig. 2). This term describes all the aftershocks triggered by mother earthquakes 
occurring within the space-time window {[£,£ + r] x S}. It is easy to show that the corresponding 
seismic rate, which can be compared with the contribution (|40p. is equal to 

MJ r ,S) = l-n + —g- II (R(r,S;x))dx, (52) 

where 



s 



(53) 

Z=l 



dz 

is the average number of aftershocks triggered within the space window S till instant r by some 
mother earthquake occurring at position x and at the time t = 0. Using relations (J26|) - (|3U|) . one 
can show that 

(R(r, S; x)) = (R(S; x)) - (R + (r, S;x)), (54) 

where (R(S;x)) is the total number of aftershocks falling inside the space domain S which are 
triggered by a earthquake source occurring at position x and at the time t = 0. (R + (t,S;x)) is 
the corresponding number of aftershocks falling with the space domain S after the instant r. 

It is easy to show that the Laplace (with respect to r) and the Fourier (with respect to x) 
transform of (R + (t,S; x)) is equal to 

(R)Ju,S;q)=n0 s (q)- ( L ) (55) 

u\\ — n(j){q) 1 — n<fi(q)<&(u) J 

Using the asymptotics (fTTj) . we can rewrite this last relation in the form 

(R) + (u, S; q) ~ (R) (S; q) ^ ~ °f f - , (56) 

T(l - 9)(cu) e + 1 - ncj)(q) 

where (R)(S; q) is the Fourier transform of the average {R{S; x)) of the total number of aftershocks 
mentioned above. Using relations (pU]l. (jSOjl . we obtain 

(R)(S;q)= ni{ f i s (q), (57) 
1 - n<p{q) 
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where Is(q) is the Fourier transform of the indicator function of the space window S. In what 
follows, we assume that S is the circular domain of radius i centered at the origin of the plane x. 
Then 

I s (q) =2tt^ J x {lq). (58) 
Eq. (JS1)|) implies that (R+(r, S;x)) is given by 

(R+(t, S; x)) = (R(S; x)) <8> H(t; x) , (59) 

where the Fourier transform of the function Ti.(r; x) is equal to 

«<-)-»(-^©> m 

Thus, the Fourier transform (with respect to x) of the sought average given by is 

1 — n<f)(q) ( t 



(R)(r,S; q)^(R)(S; q) 



I -En 



(61) 



1 — n V ci , 

Using expression ljfil)l . we construct Figure 4 which plots (R(t,S;x)) for different values of r/ci, 
in order to illustrate its convergence to the average (R(S;x)) of the total number of aftershocks 
falling inside the area S. 

As can be seen from figures 3 and 4, it follows from ()6*T|) and from the properties of Mittag-Leffler 
functions that, if the large time window approximation ()51|) holds, one may use the approximate 
equality 

(R(t,S;x))~(R(S;x)) . (62) 
In this large time window approximation, the relative rate (|S2*j) is transformed into 

n(l-n) f°° 4>{g) =2, 



As the space domain S increases in size, A/in(<S) increases towards 1. Figure 5 plots A/j n («5) using 
the space propagator <p{x) given by Q for different values of the exponent r\. 

V. LARGE TIME WINDOW APPROXIMATION 

The analysis of the previous section gives us the possibility to explore the probabilistic properties 
of the number of events in given space-time windows, in the regime where the large time window 
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approximation (|5Tj) holds. If the time duration r of the space-time window is sufficiently large, the 
previous section has shown that the statistical averages and the seismic rates become independent 
of r. It seems reasonable to conjecture that the GPF Q(z, r,S;x) of the total number of aftershocks 
triggered by some earthquake source inside the space domain S until time r coincides approximately 
with the saturated GPF Q(z, S; x) of the total number of aftershocks triggered by some earthquake 
source inside the space domain S. Within this approximation of large time windows, the effect of 
aftershocks triggered by earthquake sources occurring till the beginning t of the time window is 
negligible. Section IV 111 below will explore in more details the applicability of this conjecture. 

Within this large time window approximation, one may ignore the first term in the r.h.s. of 
Eq. (122)) and replace Q(z,t,S; x) by Q(z,S; x) in the remaining terms. As a result, Eq. (J22j) takes 
the following approximate form 

oo 

L(z, r, S) ~ r J J [I- Q(z, S; x)} [1 - I s {x)]dx + r JJ[1 - zQ(z, S; x))dx , (64) 

-oo S 

where Q(z,S; x) is the solution of Eq. with (|3Uj) or, equivalently, is the solution of 

= G[0®0-(1- z)I s e <8> 0] . (65) 
where the function G is defined in (II 7|) . 



A. Factorization procedure 

To find a reasonable approximate expression for the sought GPF Q(z,S;x), notice that if £ ^> d 
(or if n is close to 1) then the characteristic spatial scale associated with the GPF Q(z,S;x) 
becomes greater than d. Therefore, without essential error, one may replace 0®0 by in (|65p. In 
addition, we take into account the finiteness of the domain S by using the factorization procedure 
of replacing the last term of the argument of the function G in (|65|) as follows: 

I s (x)Q(z,S; x) <g> <f)(x) ~ Q(z, S; x) ps(x) , (66) 

where ps(x) remains to be specified. We will show below that ps(x) may be interpreted as the 
overall fraction of aftershocks, triggered by a mother earthquake at position x, which fall within the 
domain S. The factorization procedure amounts to replacing a convolution integral by an algebraic 
term. This factorization approximation is a crucial step of our analysis and will be justified further 
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below. As a result of its use, the nonlinear integral equation (165(1 transforms into the functional 
equation 

Q = G[(l + (z-l) Ps (x))Q]. (67) 

It is easy to show that, if relation (JBTj) holds, then the average of the number of aftershocks 
corresponding to it is equal to 

(R) = ^-p s . (68) 
1 — n 

In the next subsection, we shall clarify what is the probabilistic sense of the parameter ps- Here, 
it is sufficient to remark that one can determine it from a consistency condition: choose ps{x) such 
that the r.h.s. of Eq. (J6"%j) is equal to the true {R(S; x)). This gives 

1 — 71 ~ ~ 1 — 71 

Ps (x) = (R(S; x)) , p s (q) = I s (q)<f>(q)- ~— . (69) 

n 1 - n<p(q) 

Figure 6 plots Ps{ x ) defined by expression (JBTJj) as a function of dimensionless distance x/£. 

One can observe that, for I 3> d, the factor ps(x) approaches a rectangular function. We can 

use this observation to help determine the statistics of the number of events in a finite space-time 

window, using the approximation ps{&) — const = p for x G S. We define the parameter p as the 

space average of Ps( x ) over the window's area S: 

P - \ J J Ps{x)dx . (70) 
s 

p is thus the average over all possible spatial positions of mother earthquakes of the fraction of 
aftershocks which fall within the space-time window S. The approximation ps{x) — const = p for 
x G S allows us to simplify the last term of (jp^j) as follows: 



t II [1 - zO(z, S; x)]dx ~ tS[1 - z&(z; p)\ , (71) 
s 

where Q(z,p) is the solution of 

e(z;p)=G[(l + {z-l)p)e{z;p)), (72) 

which is derived from equation (|67|). 

Complementarily, as can be seen from figure 6, Ps{ x ) is small outside the window space domain 
S. It implies that, outside S, one may replace Eq. (JBTj) by its linearized version. As a result, we 
get 

1 - Q(z,S; x) ~ — -(1 - z)p s {x) . (73) 
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Therefore, the first term in the r.h.s. of Eq. (J54*j) transforms into 

oo 

r JJ[l-e(z,S ] x)][l-I s {x)]dx~qTS T ^(l-z) , (74) 



where 

n = 

s 

Taking into account that, due to (J69j) 



oo 



q=l II p s {x)[l-I s (x)]dx. (75) 



Ps(x)dx = S , (76) 



we obtain 

q ~ 1 — p . (77) 

Putting all these approximations together allows us to rewrite Eq. (|64j) in the form 

n 



L(z,r,S) ~ rS 



[l-p)(l-z) + l-zQ(z;p) 



(7f 



1 — n 

In what follows, we shall select a value of the parameter p which takes into account the finiteness 
of the window's spatial domain S, to better fit empirical data on the statistics of seismic rates in 
finite space-time bins. 

B. Probabilistic meaning of the factorization approximation (|66|) leading to (|67|) and 
(1531) 

The factorization approximation (jESj) has the following implication. Calling O(z) the GPF of 
the total number of aftershocks triggered over the whole space by some earthquake source, one can 
then determine the corresponding GPF Q(z, S; x) taking into account the finiteness of the space 
domain S by using the relation 

e(z,S;x) = e(q s + zp s ), q s (x) = 1 - p s (x) . (79) 

This expression f!79|) has the following interpretation. Let the above mentioned earthquake source 
triggers r aftershocks. Then, the number of those aftershocks which fall into the space domain S, 
is equal to 

R m (S;x\r)=X 1 +X 2 + --- + X rj (80) 
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where {X\, . . . , X r } are mutually independent random variables equal to 1 with probability p$ and 
with probability q$ = 1 — Ps- Thus, p$ is the fraction of the aftershocks which fall into the 
domain S. 

The corresponding expression (|79|) can be interpreted as follows. It gives the exact solution 
for the GPF of some specific space-time branching process, such that the pdf f(y; x) of the space 
positions y of each aftershock is the same for all aftershocks and depends only on the position x 
of the earthquake source. For this problem, we have 

oo 

p s (x) = If f(y;x)I s (y)dy. (81) 

— oo 

In the general case, the relation (f?9"|) offers the possibility, at least semi-quantitatively, to describe 
the characteristic features of the space-time branching processes, by using the probability ps as an 
effective independent parameter of the theory. 

Let us mention a few useful consequences of the relation (|79|). It implies that the probability 
that r aftershocks fall into the spatial domain S is equal to 

oo 

P(r,S;x) = y £P(k) B(k,r,S;x), (82) 

k=r 

where P{k) is the probability that some earthquake source triggers k aftershocks and 

B(k,r,S;x)=(*)p r s q k s - r . (83) 

This binomial probability B(k, r, S; x) is nothing but the conditional probability that, if the mother 
earthquake triggers k ^ r aftershocks then, r of them will fall into the spatial domain S. If r ^> 1, 
expression (|HH|) can be approximated by its well-known Gaussian asymptotics 

(r - kp s ) 
Zkp s q s 

If, in addition, P(k) decays slowly, for instance if it has a power asymptotic for k — > oo, then 
expressions ()82j) and ()84|) imply the asymptotic relation 



B(k,r,S;x) = exp 
y/2-Kkps qs 



4) 



P(r,S;x) ~ — P { — ) (r-Kx>). (85) 
Ps \PsJ 

In view of this last relation ([85)1 . it seems reasonable to assume that the asymptotic behavior of 
the probabilities of the number of events for r ^> 1 are the same for the case of a finite S (p$ < 1) 
and for an unbounded one (ps — 1). 

20 



C. Large space-time windows 



In order to get more insight into the properties of the statistics of seismic rates in finite space- 
time windows, it is useful to study the statistics of seismic rates in the limit where the space and 
time windows are large. In this case, A/jn in (J63|) and p in (J7U|) are both close to 1 and one may 
replace Eq. (J78j) by 

L(z,t,S)~tS[1-zB(z)], (86) 
where Q(z) is solution of the functional equation 

0(z) = G[zG(z)} . (87) 

Accordingly, the GPF of the number of events in a (large) space-time window as given by 1)31)) 
takes the form 

e sp (z,p) = f( P [i-ze(z)}). (88) 

Here and everywhere below, 

p=(g)rS. (89) 

Knowing the GPF Q sp (z,r;S), the probability P sp (r; p) of event numbers r is obtained from 
the formula 

P, p(r;p) i */(,n -*>(»)]) . (90) 



r! dz r 
Equivalently, the integral representation of (J90|) reads 



z=0 



^s P (r; p) = ^-f f(p [1 " ze{z)))^ , (91) 
c 

where C is a sufficiently small contour in the complex plane z around the origin z — 0. 

The main difficulty in calculating P sp (r; p) comes from the fact that the GPF O(z) is defined 
only implicitly by Eq. (|87|). To overcome this difficulty, we rewrite the integral in ()91|) in the 
following equivalent form 

P, p(r;p) = J_/ ^H :; e W ]) (r>0) (92) 

C 

and use the new integration variable y = zQ(z). Expression ()87|) shows that the inverse function 
of y is z — y/G{y). As a result, the equation (}9*2j) transforms into 

Psp(r;p)= 2^/ Gr(2/)g(2/;/?) l' (93) 
c y 
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where 



(94) 



p dz 

and C is a contour enveloping the origin y = in the complex plane y. One may interpret Q(z; p) 
((SH) as the GPF of some random integer R p such that (i? p ) = p. 

Notice that Eq. has a simple probabilistic interpretation. Indeed, it follows from that 



P sp (r; p) = ^ Pr {i? p + R{r) = r - 1} , 
where i? p is a random integer with GPF Q(z; p) given by (J9~2j) while 



R(r) = R 1 + R 2 + --- + R r 



(95) 



(96) 



where {R\, R2, • • • , -Rr, • • • } are mutually independent random integers with GPF G(z). This im- 
plies that the probabilities of each such random variable Ri, i = 1, r, has the power law asymp- 
totics (fPQl. 

This remark provides a simple analysis of the asymptotic behavior of the probabilities P S p(^; p) 
for r 3> 1, by using expression f)95j) . For 1 < 7 < 2, the asymptotics of the probability P(A;; r) that 
the sum (|9l)|) is equal to k goes to, for large r, 

A; — nr \ 



P(Jfe|r) 



er 



er 



1 1/7 



where 



n- 



7 , 



r(i- 7 ) 



and V ; 7( a; ) is the stable Levy distribution with the two-sides Laplace transform 

/oo s 
ips{x)e~ sx dx = e s . 
-00 

It is known that 



x 



-S-l 



T(-5) 



(X — > OO) , 



MO) 



1 



*r(i-i: 



(97) 



(9* 



(99) 



(100) 



One can calculate V^^) for any 1 < 5 < 2, using, for instance, the following integral representation 



1 f°° 

il) S (x) = - exp 

7T JO 



— W + UX COS 



7T 



sin 



ux sin 



7T 



(101) 
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For some numerical illustrations, we will use the case 5 = 3/2 for which the following analytical 
expression is available 



^) = ^ 5 [r(l) 1 Fi(i,l,^)-^r(|) 1 F 1 ( 



'7 4 4^ 
6' 3' 27 



(102) 



For r>/), one can neglect the random integer R p in the r.h.s. of Eq. ()95|) and obtain from (|95j) 
and ()97j) the following asymptotic formula 



p 



If 1 — n <C 1 (the branching is close to but not exactly critical), Eq. ()103j) predicts the existence 
of two characteristic power laws in the dependence of the probabilities P sp (r;p) with r, a result 
already derived in 



'1 — n)r — 1 



er 



,1/7 



(r > p) . 



(103) 



id- 



1. For 



then, 



2. For 



1 r r* , with r* 



1 v 7/(7-1) 

1 \ C V(7-1) 



P S p{r;p) ~ r 



1 — n 



-1-1/7 



we recover the original power law (|19J1 of the number of first generation aftershocks 



P sp (r;p) ~ r 



-1-7 



(104) 



(105) 



(106) 



(107) 



For values of the parameters 7 = 6/0 and n which are typical of real seismicity modeled by 
aftershock triggering processes, the cross-over number r* separating the two power laws can be 
very large. For instance, for 7 = 1.25 and n = 0.9, we obtain r* ~ 10 4 . 



D. Prediction of the distribution of event numbers for large time windows 

Starting from the general expression (j3~T|) of the GPF sp (z, r; S) with the approximation (f7S|) 
for L(z, t, S) and using the relationship between the probability P sp (r; p,p) and its GPF similar to 
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expression (j9~Uj) and its integral representation similar to ([91)1 . we obtain the following expression 
valid in the limit of large time windows 



P sp (r;p,p) = —ff[p 
c 



n 



n 



-{l-p){l- z) + l- ze{z;p) 



dz 



yr+l 



(108) 



where Q(z;p) is the solution of the functional equation (|72*|) . Similarly to the change of variable 
used to go from (|9*2|) to (J9*H|) . we introduce the new integration variable 



y = (l + (z-l)p)Q{z;p) 



z{y) = - 



p 



(109) 



V \G(y) 

By construction of y, Q(z;p) = G(y) which allows us to obtain the following explicit expression 

P sp {r;p,p) = ^-x 

dZ{y) dy 



f P 



n 



. 1 — n 



l-p)(l-Z(y)) + l-Z(y)G(y) 



(110) 



dy Zr+i(y)' 



This expression (111 U j) allows us to make a precise quantitative prediction for the dependence of 
the distribution P sp (r;p,p) of the number r of earthquakes per space-time window as a function 
of r, once the model parameters n, j,p, S and p are given. We note that Pisarenko and Golubeva 
lave shown that the distribution of numbers has the same tail as the distribution of seismic rates 
^]. Thus, for the tails, our results can be interpreted either as statements on the distribution of 
realized numbers of earthquakes or on the distribution of average seismic rates. 

We now turn to a brief description of the data analysis and of their fits with (jllOj) . 



VI. EMPIRICAL ANALYSIS AND COMPARISON WITH THEORY 

We use the Southern Californian earthquakes catalog with revised magnitudes (available from 
the Southern California Earthquake Center) as it is among the best one in terms of quality and 
time span. Magnitudes Ml are given with a resolution of 0.1 from 1932 to 2003, in a region from 
approximately 32° to 37°N in latitude, and from —114° to —122° in longitude. In order to maximize 
the size and quality of the data used for the analysis (to improve the statistical significance), we 
consider the sub-catalog spanning the time interval 1994 — 2003 for Ml > 1.5, which contains a 
total of 86, 228 earthquakes. The completeness of this sub-catalog has been verified in the standard 
way in j3] by computing the complementary cumulative magnitude distribution for each year from 
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1994 to 2003 included. The stability of the linear relationship of the logarithm of the number as a 
function of magnitude Ml for Ml > 1.5 is taken as a diagnostic of completeness. 

The spatial domain are covered by square boxes of (L = 5 km) x (L = 5 km), which gives us 
16046 spatial bins, with many of them being actually empty. This size is a compromise between 
conflicting requirements. On the one hand, a smaller spatial resolution should not be used due 
to the errors of localization of earthquake epicenters, which are of this order. On the other hand, 
increasing L decreases very fast (~ l/-^ 2 ) the total number of boxes with which the distribution of 
the number of events can be constructed. Four different sizes for the time window are considered: 
dt = 1 day (3652 temporal bins), dt = 10 days, dt = 100 days and dt = 1000 days. Combining 
the space and time windows leads to space-time windows or bins of size L 2 x dt. For instance, 
for dt = 1 day, there is a total of 54669 spatio-temporal bins with at least one event, with 4298 
non-empty spatial bins. 

Figures 7-10 plot the empirical probability density functions Pdata(^) of the number r of earth- 
quakes in the space-time bins described above. The straight line in Figs. 8-10 is the best fit with 
a pure power law 

Pdata(r) ~ l/r 1+ " (111) 

over the range 10 < r < 100. The estimation for /x is found stable across different time windows, 
since the fitted values are /x = 1.65 for dt = 100 days, /x = 1.75 for dt = 10 days and /x = 1.60 
for dt = 1 days. However, one can also observe at the same time that the pdf becomes more and 
more curved in the larger portion of the bulk as the size dt of the time window is increased. This 
behavior can be explained by our theory as we know describe. 

First, all parameters are let free to adjust optimally. The parameters are the branching ratio n 
defined by (JBJ), the exponent 7 defined in (J3J), the exponent S defined in ()15|). the parameter p$ ~ p 
coming from the factorization procedure (j66J) and given by (JBTJj) and ()70|) (which is roughly equal to 
the overall fraction of aftershocks triggered by sources in the domain S which fall within the same 
domain S), and the average number p of spontaneous earthquake source per space-time bin defined 
in (|89)1. The theoretical curves shown in Figures 7-10 are obtained by a numerical integration of 
()110|) for the set of parameters n = 0.96, 7 = 1.1, p = 0.25, 5 = 0.15 and p = 0.0019 dt days 
with dt in units of days (thus equal to 1000 for Fig. 7). Note that a given p for a given space-time 
window [t, t + t] x S translates into the following average number of events inside that window: 
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(R sp (t,S)) ~ ~ 0.02 — 0.16 x dt (in units of days), taking the mean value n = 0.96. These 
parameters are chosen to best fit the empirical pdf for the largest time window dt = 1000 days 
shown in Fig. 7. They are then frozen and the theoretical distribution is recalculated with formula 
(|110|) with no adjustable parameters for the three other cases dt = 100 days (Fig. 8), dt = 10 days 
(Fig. 9) and dt = 1 day (Fig. 10). The theory is thus able to account simultaneously for all the 
considered time windows, with no adjustable parameters for the three smallest time windows. 

Second, we test for the sensibility of the parameters. We actually obtain practically the same 
quality of fits for all four values dt = 1, 1000 days by varying 5 and 7 in the ranges 0.1 < 5 < 0.2 
and 1 < 7 < 1.5, by properly adjusting the other parameters adaptively. We do not show the 
corresponding theoretical curves as they are essentially equally good to fit the data within the 
empirical noise. The choice 5 close to zero is consistent with the choice of the Cauchy distribution 
as a proxy for the heterogeneity of the spatial distribution of spontaneous earthquake sources 



and empirical analysis 



inferred for the stress field and deduced from previous theoretical 

n 

of earthquake sources 24J. The strong sensitivity of the fits with respect to the fractal structure 
of the spontaneous sources quantified by the exponent 5 is a surprising but positive return of this 
work. We did not expect a priori that the distribution of seismic rates would teach so much about 
the heterogeneity of the seismic active regions. But there is a lack of sensibility with respect 7. 
Thus, the distribution of seismic rates cannot be used to constrain the productivity parameter a 
thr oug h 7 = b/a) other than by confirming the range of previous estimations of its value: a ~ 1 

aJ0i, 0.5 < a < 1 0Q,Q- 

The sensitivity of the fits with respect to the branching ratio 
n and to the overall fraction p of aftershocks which fall within the domain S requires a special 
discussion which will be reported elsewhere. 

We have also used functional forms for fs{x) other than ()15|) to describe the pre-existing hetero- 
geneity of spontaneous earthquake sources, such as half-Gaussian, exponential as well as different 
variants of power laws. Overall, we find that we need fs{x) to have a power law tail close to the 
Cauchy distribution in order to get a reasonable fit for all four time windows. This shows that the 
ETAS model as well as any other model of this class of triggered seismicity need to be general- 
ized to account for a pre-existing heterogeneity of the crust, which controls the occurrence of the 
spontaneous earthquake sources. 

We would also like to stress that, according to our theory, the value of the exponent \i ~ 1.6 
used in (jlllj) to fit the tails of the distributions shown in Figs. 8-10 is describing a cross-over rather 
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than a genuine asymptotic tail. Recall that the distribution of the total number of aftershocks has 
two power law regimes ~ l/r 1+ ~ for r < r* ~ 1/(1 — n) 7 ^ 7-1 ) and ~ l/r 1+7 for r > r„ jis| . 
The existence of this cross-over together with the concave shape of the distribution at small and 
intermediate values of r combine to create an effective power law with an apparent exponent /i ~ 1.6 
larger than the largest asymptotic exponent 7. We have verified this to be the case in synthetically 
generated distributions with genuine asymptotics exponent 7 = 1.25 for instance, which could be 
well fitted by // ~ 1.6 over several decades. We note also that Pisarenko and Golubeva with 
a different approach applied to much larger spatial box sizes in California, Japan and Pamir-Tien 
Shan, have reported an exponent fi < 1 which could perhaps be associated with the intermediate 
asymptotics characterized by the exponent I/7 < 1, found in our analysis By using data 

collapse with varying spatial box sizes on a California catalog, Corral finds that the distribution 
of seismic rates exhibits a double power-law behavior with fi « for small rates and /i ~ 1.2 
for large rates [2]. The first regime might be associated with the non universal bulk par of the 
distribution found in our analysis. The second regime is perhaps compatible with the prediction 
for the asymptotic exponent /i = 7. 



VII. THEORETICAL TESTS OF THE THEORY USING STATISTICS CONDITIONED 
ON GENERATION NUMBER 

In our theoretical development to obtain the prediction ()110|) that could be compared with 
empirical data, we have been obliged to make two approximations: (1) assuming that the duration 
r of the time window [t, t + r] is sufficiently large (i.e., the inequality (JoTj) holds), we have replaced 
the GPF Q(z, t,S;x) by its asymptotics 0(z, S; x); (2) we have used a factorization procedure to 
take into account quantitatively the finiteness of the spatial window S. 

In this section, we attempt to clarify further the domain of application of these two approx- 
imations by testing them on other event statistics conditioned on fixed number of generations. 
Numerical calculations of the exact PDF are compared with the approximations. 



27 



A. Large spatial windows 



Let us consider the statistics of aftershocks triggered over the whole space during the first k 
generations by some mother event. The corresponding GPF of the number of aftershocks triggered 
in the course of k generations is defined by the following iterative recurrence equation 

Q k (z) = G[zQ k ^(z)] , ei(*) = GO) . (112) 

One can calculate the corresponding probabilities of aftershock numbers by using the Cauchy 
integral 

P>(r) = ±fe i {z)£ [ . (113) 
c 

Furthermore, we can make use of the knowledge that, as k — > oo, the GPF Qk(z) converges to the 
asymptotic GPF Q(z) which is the solution of Eq. (|S7jl. It is easy to show that one can calculate 
the corresponding probabilities using an equality analogous to 

P(r) = 7^ r / G r+1 (y)-^- . (114) 

Figure 11 shows the distribution Pfe(r) of aftershocks numbers, obtained by a numerical calculation 
of the integral ()113j) for different values k = 1,2, 3, 5, 8 and for 7 = 1.25 and in the critical 
case n — 1. It also shows the corresponding asymptotic distribution for k —>■ +00 obtained by 
integration of (|114|) . Note that, even for in this critical case, the distribution for k = 8 generations 
is already almost undistinguishable from the asymptotic distribution including an infinite number 
of generations, at least for r ^ 250. Figure 12 clarifies further the convergence rate by plotting the 
ratio 

P k (r) = ^ (115) 

for different values k. 

The information on the number of generations necessary to reach the asymptotic regime gives 
us the possibility of estimating the corresponding characteristic time beyond which the asymptotic 
distribution P(r) becomes a good approximation of Pfc(r). Let T(k) denote the random time at 
which a k-th generation aftershock is triggered. It is equal to 

T{k) = n + r 2 + • • • + r k , (116) 
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where {n, T2, ■ ■ .Tk, ■ ■ ■ } are mutually independent random waiting times between the occurrence 
of a mother earthquake and one of its first-generation aftershock. We define the u-th waiting time 
quantile t(u, k) of generation k by 

Q[t(u, k), k] = Pr{T(k) > t(oj, k)} = to . (117) 

Thus, 1 — u is the probability that the duration of any chain of k successive generations of triggered 
aftershocks is smaller than t(u, k). Choosing some confidence level (for example 1 — u — 0.9), one 
may assert that, during the time t(u>, k), all aftershocks of the k-th generation have already been 
triggered. Let k — h* be such that the corresponding probability PkS r ) ls close to the asymptotic 
P(r). Then, one may interpret 

U = t(u,k # ) (118) 

as an estimation of the characteristic time for the validity of the asymptotic distribution P(r). 

The asymptotic expression for the probability Q(t, k) defined by (|117|) for k ^> 1 can be deter- 
mined by using the fact that the terms of the sum (|116|) are determined by Omori's law (jHJ) with 
< 6 < 1. For k ^> 1, Q(t, k) is asymptotically close to 

fl(t -* ) = F ' Ur ( i -»)]■/» ) ' (119) 

where 

F e (x) = / Mv)dy (120) 



and <po(x) is the one-sided Levy stable distribution defined by the Laplace transform 

roc „ 

<p e (u)= V e(x)e~ ux dx = e~ u . (121) 







In particular 

F 1/2 (x) = erf (-L=) . (122) 



2^J ' 



The following asymptotic behavior holds 



^fr)- r(i-g) ( x>>1 )- ( 123 ) 

Substituting (|119j) and (|123|) into ()117jl . we obtain the following estimation for t* defined by (|118|) 

(k \ 1/e 

^d- • (124) 
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For instance, for 9 = 1/2, k* = 8, to = 0.1 and c = 2 minutes, then £* ~ 9 days. Note that £* is 
highly sensitive to the value of 9. Indeed, = 1/3 (resp. 2/3) with all other parameters being the 
same gives ~ 700 days (resp. £* ~ 1 day). 

Expression (|124j) is related to condition (J51)) (and actually improves on it) as follows. The 
condition ()51|) with (joTIj) can be expressed by introducing, similarly to the reasoning leading to 
f!124|) . some small threshold to <C 1 such that (JoTj) translates into A^ou^t) % Correspondingly, 
we can introduce r* such that, if r > r*, then the duration of the window is large at the u level: 

«=^teV- (125) 



r(l-0) Vr 
Using (US)), we get 

/ \ 1/0 

r*^c — ^- . (126) 
Vo;(l-n)y 

The two expressions ()12fi)) and ()124j) have a similar structure. The only difference is that the 
characteristic generation number /c* is replaced by the factor n/(l — n). For n not too close to 1, 
n/(l — n) gives a not unreasonable estimation of k*. For n close to 1, expression ()124|) should be 
preferred as it provides an improvement to ()51|) based on the calculation of quantiles rather than 
on the mean rate behavior. 



B. Finite spatial windows 

A natural generalization of the iterative procedure ()112j) for finite spatial window S allows 
to estimate the corresponding aftershock statistics. Consider an earthquake occurring at point x. 
Then, the pdf of the space positions y of an aftershock of the k-th generation is given by <fik{y — x), 
where 

Mv) = <t>(y) <t>(y) (127) 

k 

is the fc-times convolution of the space propagator <f)(y) (one example is given by ©)• Correspond- 
ingly, the probability for an aftershock of the k-th generation to fall into the space window S is 
equal to 

p k (S;x) = J J <f) k (y - x)dy . (128) 
s 
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Provided these probabilities are known, one can determine the GPF Qk{z,S; x) of the number of 
aftershocks of the k generation occurring in the spatial domain S by using the following iteration 

Q 1 (z,S;x)=G( Pl (z-l) + l) , 
2 (z, S;x) = G [( Pl (z - 1) + 1) G (p 2 {z - 1) + 1)] , (129) 
e 2 (z, S;x) = G [( Pl (z - 1) + 1) G [(p 2 (z -1) + 1)G i(p 3 (z - 1) + 1)]]] , 

and so on up to the order k. Then, the distribution of the number of aftershocks of the k generation 
is given by 

I f dz 
P k (r, S\x) = —j Q k (z, S;x)—. (130) 

c 

These expressions are general and hold for any space propagator 4>(y). Let us now specialize to 
the form (jUJ) for the spatial propagator (p(y), with rj = 1, 

kd 

MX) = 2n(x 2 + ¥dW* (131) 

corresponding to an asymptotic l/|x| 3 decay law often argued on the basis of the shape of the 
elastic Green function in a three dimensional space. From (|128jl . we then have 

= (132) 

where 

^/ x 2 r u „ / Avs \ sds ,,„. 

V{u,v) = - El- -, . 133 

7T Jo \1 + (v + s) 2 J [l + ( v - s) 2 ]^fl + {v + s) 2 

Here, E(m) is the complete elliptic integral 

«r/2 / 

E(m)= yi - msin 2 erfe. (134) 
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Figure 13 shows the distributions ()132|) for £ = lOd (recall that £ is the radius of the assumed 
circular domain S centered on the origin, which has been used in (|58p) and for different positions 
x of the mother earthquake, given by x/£ = 0; 0.4; 0.6; 0.8; 1; 1.2; 1.4; 1.6 from top to bottom. 
The separation by the curve for x/£ — 1 into two families has a simple explanation. For x/£ < 1, 
the mother event lies within the spatial domain S of interest and it is thus counted as generation 
0. Its immediate aftershocks are most probably adjacent to it and thus have a large probability 
to also fall within S. As the number k of generation increases, aftershocks diffuse away and are 
less and less likely to fall within S. In contrast, for x/£ > 1, the mother earthquake falls outside 
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S. Therefore, there is not event at the zeroth generation in S, hence the curves start from zero. 
The first generations of aftershocks which are most likely to be nearby the mother earthquake fall 
rarely within S. Only as aftershocks of higher generation levels develop and diffuse away from the 
mother earthquake, can they invade S. Of course, at large generation numbers, the aftershocks 
diffuse away from any finite spatial domain, explaining the decay of Pk{<S', x) to zero for large fc's. 

Figure 14 plots the asymptotic distribution P(r,S; function of the number r of after- 

shocks for 7 = 1.25, n = 0.99, for different values of the radius £ of the disk S. The mother 
earthquake is assumed to occur at the origin, that is, at the center of the disk S. P(r,S;x) is 
obtained by using ()130)) for k = 25 generations, which is certainly a very good approximation to 
P(r,S;x) = Pk^ +00 (r,S;x). Figure 15 plots P(r, S;x) as a function of the number r of after- 
shocks, at fixed £/d = 10 for various positions of the mother earthquake, for the same parameters 
7 = 1.25, n = 0.99. One can observe that, when the mother earthquake is inside the disk S 
(x = 0.8£; 0.6£; 0A£; 0.2£; 0), the corresponding distributions are close to each other as predicted 
in Section IV Al When the mother is outside the disk S (x = 1A£; 1.2£; 1), the distributions differ 
from the previous case and are significantly smaller. This gives additional support in favor of 
the linear approximation (fT3*j) . which we used in Section fV Al More precisely, these properties 
result directly from the analysis of section IV At in which we notice below equation ()69|) that, for 
£ ^> d, the factor ps{x) approaches a rectangular function, which leads to the approximation 
Ps{x) ~ const = p for x £ <S. This leads to the natural assumption that the GPF Q(z,S; x) is 
almost the same for all interior source positions x 6 S. This means in turn that the corresponding 
distribution P(r,S;x) should be the same for all x £ S. This remarkable fact is illustrated in 
figure 15 in which the curves for the interior source probabilities merge. 

C. Testing the factorization approach 

We can now test the factorization approximation developed in section IV Al to take into account 
the finiteness of the space window S by comparing it with the approach in term of the statistics 
over successive generations of the previous section. We thus compare the asymptotic distribution 
P(r, S, x — 0), obtained by calculating the integral ()130|) for a large enough generation number k 
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(k = 25 is found to be sufficient), with the factorization approximation 

f j dGjy) G r {y)dy 

p ^-2^f^r [y -(i- P )G( V )v (135) 

with an appropriate value of the parameter p = Ps, defined as the fraction of the aftershocks 
which fall into the domain S. Figure 16 shows this comparison for four different values of E/d = 
10; 7.5; 5; 2.5. The upper curve in each panel is the distribution (|114j) for an infinite domain 
£/d — +oo, as a reference. There is some discrepancy between P(r,S,x = 0) and P(r,p) given 
by (|135p . The main difference is that the true distribution P(r,S,x = 0) decays faster than the 
factorization approximation P(r,p) for large r. We recall that, due to ([85)1 . the asymptotic behavior 
of the distribution P(r,p) obtained under the factorization approximation is the same as for an 
infinite domain given by (J114)) . Nevertheless, Figure 16 shows that an appropriate choice of the 
parameter p = ps allows us, at least semi- quantitatively, to take into account the finiteness of the 
area S. 



VIII. DISCUSSION 



We have presented a general formulation in terms of generating functions of the space-time orga- 
nization of earthquake sequences, in the framework of general branching processes. We have applied 
this approach to the ETAS (Epidemic- Type Aftershock Sequence) model of triggered seismicity. 
In view of the formidable difficulty in obtained exact solutions to the nonlinear integral equations 
involving the generating functions, we have developed several approximation schemes which have 
been tested by comparison with exact numerical calculations. We have used the corresponding 
predictions to fit the distribution of seismic rates in four finite space-time windows in a California 
seismic catalog. The space-time windows differ by their time interval going from dt = 1 day to 
dt = 1000 days. The fits have been found to account satisfactorily for the empirical observation. 
In particular, we have adjusted the parameters of the theory on the largest time window dt = 1000 
days and have then used these frozen parameter values in the theory to calculate the distribution 
for the smaller time windows. This tests the rigidity of the theory to account simultaneously for 
the distributions at different time scales. In this process, we have found it necessary to augment the 
ETAS model by taking account of the pre-existing frozen heterogeneity of spontaneous earthquake 
sources. We have discussed the physical justification of this generalization in terms of pre-existing 
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stress and fault networks, which constrain the form of the pre-existing heterogeneity. Our find- 
ings have also important implications to assess the quality of models developed to forecast future 
seismicity, and suggest to re-examine current procedures which assume Poisson statistics in the 
construction of likelihood scores. 
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P 1 (r) 




Fig. 1: Plot of the probabilities ()18|) and their power law asymptotics ()19(l for the infinite 
variance case 7 = 1.25 and for a finite variance case 7 = 3. 
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Fig. 3: Plots of the exact rate A/out CO) its fractional approximation (|46|) (which actually 
coincides with the exact value) and its asymptotic approximation (|50ft obtained from ()48|) 
(dashed line), for n = 0.9 and 6 = 1/2. 
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Fig. 4: Plots of {R(t,S;x)) for 7/ = 1, £ = 10 d, 9 = 1/2 and for r = ci; 5ci; 20 ci. Dashed 
line is the plot of average of total number of aftershocks, falling into the circle S. 
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Fig. 5: Plots of the relative rate A/i n (5) given by ()63|) as a function of the dimensionless size 
£/d of the space domain, for different exponents r] = 1; 2; 3 of the space propagator. 
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Fig. 6: Plots of the self-consistent factor ps(x) defined by (|69j) for £ = 10 d and for n 
0.5; 0.8; 0.9. 
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Fig. 7: Empirical probability density functions -Pdata( r ) of the number r of earthquakes in 
the space-time bins of size 5x5 km 2 and dt = 1000 days. The continuous line is the fit of 
this data with formula (|110|) for the set of parameters n = 0.96, 7 = 1.1, ps = 0.25, 5 = 0.15 
and p = 0.0019 km xlOOO days. 
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Fig. 8: Empirical probability density functions -Pdata( r ) of the number r of earthquakes in 
the space-time bins of size 5x5 km 2 and dt = 100 days. The straight line is the best fit with 
a pure power law -Pdata( r ) ~ l/r 1+At over the range 10 < r < 100, which gives \i = 1.65. The 
continuous line is the theoretical prediction using formula (jllOj) with dt = 100 days and with 
no adjustable parameters, as the parameters are fixed to the values adjusted from the fit in 
Fig.7. 
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Fig. 9: Empirical probability density functions -Pdata( r ) of the number r of earthquakes in 
the space-time bins of size 5x5 km 2 and dt = 10 days. The straight line is the best fit with 
a pure power law Pdata( r ) ~ l/r 1+At over the range 10 < r < 100, which gives fi = 1.75. The 
continuous line is the theoretical prediction using formula (|110|) with dt = 10 days and with 
no adjustable parameters, as the parameters are fixed to the values adjusted from the fit in 
Fig.7. 
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Fig. 10: Empirical probability density functions -Pdata( r ) °f t ne number r of earthquakes 
in the space-time bins of size 5x5 km 2 and dt = 1 day. The straight line is the best fit with 
a pure power law -Pdata( r ) ~ l/r 1+At over the range 10 < r < 100, which gives \i = 1.60. The 
continuous line is the theoretical prediction using formula (j!10|) with dt = 1 day and with 
no adjustable parameters, as the parameters are fixed to the values adjusted from the fit in 
Fig.7. 
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Fig. 11: Distributions -Pfc(r) given by for 7 = 1.25, n = 1, for different numbers of 

generations k = 1; 2; 3; 5; 8. The upper curve corresponds to asymptotic distribution P(r) 
given (|114|) corresponding to k — > +00, while lower one corresponds to the distribution Pi(r) 
given by (|18|) of the number of aftershocks of first generation. 
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Fig. 12: Ratios as a function of event numbers r for different values of the genera- 

tion number k, demonstrating the convergence of the distributions Pfe(r) to the asymptotic 
distribution P(r). Bottom to top: k = 2; 4; 6; 8; 10. 
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Fig. 13: Plots of the probabilities P k (r, S; x) given by gSSj) for 7 = 1.25, n = 0.99, £ = lOd 
and for different positions of the mother earthquake: x/l = 0; 0.4; 0.6; 0.8; 1; 1.2; 1.4; 1.6. 
Recall that £ is the radius of the assumed circular domain S centered on the origin. The two 
families of curves separated by the central one for xjl = 1 are explained in the text. 
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Fig. 14: Plots of the distribution P(r,S;x) of the total number r of aftershocks falling 
within the disk S for a mother earthquake at the origin x = and for different values of the 
circle radius I. Bottom to top: l/d = 1; 3; 5; 10; 20. 
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Fig. 15: Plots of the distribution P(r, S; x) of the total number r of aftershocks falling within 
the disk S for different positions x of the mother earthquake at fixed disk radius £/d = 10. 
Bottom to top: xjl = 1.4; 1.2; 1; 0.8; 0.6; 0.4; 0.2; 0. The upper curve thus corresponds to 
an infinite disk t = +oo. 
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Fig. 16: Comparison of the asymptotic distribution P(r,S,x = 0), obtained by calculating 
the integral IJIHUJI for a large enough generation number k (k = 25 is found to be sufficient), 
with the factorization approximation P(r,p) given by (jl35|) . where p = ps is defined as the 
fraction of the aftershocks which fall into the domain S. The upper curve in each panel is 
the distribution (|114[) for an infinite domain ljd= +oo, as a reference. The different panels 
correspond to t/d = 10; 7.5; 5; 2.5, with 7 = 1.25, n = 0.99, x = 0. 
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